Euler
ExamplesThe Newton algorithm replaces the function with its tangent, and takes the zero of the tangent as an approximation of the zero of the function. What if we use the Taylor polynomial of degree 2?
Let us start with the ordinary Newton Algorithm. First we compute the tangent in x=a of a general function f(x).
>function p(x,a) &&= f(a)+'diff(f(a),a)*(x-a)
d
(-- (f(a))) (x - a) + f(a)
da
Find the zero of the tangent. This is the Newton operator, which maps a to na for a better approximation of the zero of f.
>&solve(p(x,a)=0,x); newton &= rhs(%[1])|expand
f(a)
a - ---------
d
-- (f(a))
da
The operator has each zero of f as a fixed point. The contraction of such an operator depends on the first derivative of the operator in the fixed point.
Th Newton algorithm is fast, since the derivative of the operator in the zero is 0.
>dnewton &= diff(newton,a)
2
d
f(a) (--- (f(a)))
2
da
-----------------
d 2
(-- (f(a)))
da
Indeed. With Maxima's atvalue preset, we can teach Maxima to compute the derivative of the operator in the zeros of f.
>&atvalue(f(x),x=a0,0); &at(dnewton,a=a0)
0
We get the following Taylor expansion for the operator.
>&taylor(newton,a,a0,2)
!
2 !
d ! 2
(--- (f(a))! ) (a - a0)
2 !
da !
!a = a0
a0 + -----------------------------
!
d !
2 (-- (f(a))! )
da !
!a = a0
It starts with (a-a0)^3, which tells us that (na-a0) is less equal than a constant times (a-a0)^2. This is called quadratic convergence.
For an example, we solve a^3=8.
>&subst(f(a)=a^3-8,newton); na &= % | nouns
3
a - 8
a - ------
2
3 a
Here is the exansion of the operator, showing quadratic convergence.
>&taylor(na,a,2,3)
3 2
(a - 2) (a - 2)
- -------- + -------- + 2
3 2
Test with 100 digits precision.
>fpprec &= 100;
In each step, we double the number of valid places.
>&at(na,a=2.2b0), &at(na,a=%), &at(na,a=%), &at(na,a=%), &at(na,a=%)
2.017630853994490358126721763085399449035812672176308539944903\
581267217630853994490358126721763085399b0
2.000153616593130364027516166650181729183342359874352201195091\
532576002725023944107721337345875938096b0
2.000000011797820607418091417206088019779634849405539582686471\
063401306518202668409505552656957631673b0
2.000000000000000069594284995035565044396576173429581450020821\
753424329143996927312275962422339266909b0
2.000000000000000000000000000000002421682251985116086491265907\
737608871913078530078031916652981710816b0
Let us try to improve this using the Taylor polynomial of degree 2.
>function p(x,a) &&= f(a)+'diff(f(a),a)*(x-a)+'diff(f(a),a,2)*(x-a)^2/2
2
d 2
(--- (f(a))) (x - a)
2
da d
--------------------- + (-- (f(a))) (x - a) + f(a)
2 da
First of all, these polynomials have two zeros. Let us take the second one. We get a somewhat ugly operator.
>&solve(p(x,a)=0,x); newton2 &= rhs(%[2])|expand
2
d 2 d
sqrt((-- (f(a))) - 2 f(a) (--- (f(a)))) d
da 2 -- (f(a))
da da
---------------------------------------- - ---------- + a
2 2
d d
--- (f(a)) --- (f(a))
2 2
da da
Since we take the larger zero of the parabola, we assume that the derivative of f at a0 is positive. We can teach that to Maxima.
>&assume(at('diff(f(a),a),a=a0)>0);
With this information, Maxima can compute the derivative of the operator near a0.
>&diff(newton2,a); &at(%,a=a0)
0
This time the second derivative is 0 too.
>&diff(newton2,a,2); &at(%,a=a0)
0
So we get cubic convergence if we compute the Taylor expansion of the operator near a0.
>&taylor(newton2,a,a0,3)
!
3 !
d ! 3
(--- (f(a))! ) (a - a0)
3 !
da !
!a = a0
a0 - -----------------------------
!
d !
6 (-- (f(a))! )
da !
!a = a0
Let us test our example.
>&subst(f(a)=a^3-8,newton2); na2 &= % | nouns
4 3
sqrt(9 a - 12 a (a - 8)) a
-------------------------- + -
6 a 2
The taylor expansion in the zero 2 of f.
>&taylor(na2,a,2,3)
3
(a - 2)
2 - --------
12
In each step, we triple the number of valid digits.
>&at(na2,a=2.2b0), &at(na2,a=%), &at(na2,a=%)
1.99932634721099925522045624728277395555757289734433479625388\
9746839036446657406931706578324785885588b0
2.00000000002547575946075408666196905903116295761445444758303\
9561889029497957818231105556768840585764b0
1.99999999999999999999999999999999862215560787428855724251039\
2199433240948356638428015183239783315089b0
However, the procedure is not worth the effort, since two ordinary Newton steps are better than one step of this algorithm, and much easier to carry out, since we do not have to compute a square root.
>na3 &= at(na,a=na)|ratsimp, &at(na3,a=2.2b0), &at(na3,a=%), &at(na3,a=%),
9 6 3
4 a + 102 a + 192 a + 256
----------------------------
8 5 2
9 a + 72 a + 144 a
2.00015361659313036402751616665018172918334235987435220119509\
1532576002725023944107721337345875938096b0
2.00000000000000006959428499503556504439657617342958145002082\
1753424329143996927312275962422339266909b0
2.00000000000000000000000000000000000000000000000000000000000\
00000029322724647898516428157063672215b0